Theoretical estimates for the largest Lyapunov exponent of many-particle systems 
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The largest Lyapunov exponent of an ergodic Hamiltonian system is the rate of exponential growth 
of the norm of a typical vector in the tangent space. For an A'^-particle Hamiltonian system, with a 
smooth Hamiltonian of the type p^-\-V{q), the evolution of tangent vectors is governed by the Hessian 
matrix V of the potential. Ergodicity implies that the Lyapunov exponent is independent of initial 
conditions on the energy shell, which can then be chosen randomly according to the microcanonical 
distribution. In this way a stochastic process V(t) is defined, and the evolution equation for tangent 
vectors can now be seen as a stochastic differential equation. An equation for the evolution of the 
average squared norm of a tangent vector can be obtained using the standard theory in which the 
average propagator is written as a cummulant expansion. We show that if cummulants higher than 
the second one are discarded, the Lyapunov exponent can be obtained by diagonalizing a small- 
dimension matrix, which, in some cases, can be as small as 3x3. In all cases the matrix elements of 
the propagator are expressed in terms of correlation functions of the stochastic process. We discuss 
the connection between our approach and an alternative theory, the so-called geometric method. 

PACS: 05.45.-|-b; 05.20.-y; 02.50.Ey 



I. INTRODUCTION 

The largest Lyapunov exponent measures the sensitiv- 
ity to initial conditions in a dynamical system. In low- 
dimensional models, the Lyapunov exponent sets a limit 
to the prediction of the time evolution of a given state 
of the system. In the high dimensional systems encoun- 
tered in thermodynamics, one abandons from the start 
a description in terms of single (microscopic) states, and 
resorts to a statistical approach. In these cases a positive 
Lyapunov exponent is usually welcome, as it is a neces- 
sary condition for the validity of the Boltzmann-Gibbs 
scenario. A Lyapunov exponent that becomes null in the 
thermodynamic limit is a signal of anomalous behaviors. 
For instance, metastable phases of some long-range inter- 
acting systems, where the Lyapunov exponent vanishes 
in the large N limit, exhibit breakdown of ergodicity, 
anomalous diffusion, and non-Maxwell velocity distribu- 
tions 1^. An extension of standard statistical mechanics 
is required for the theoretical explanation of such phe- 
nomena [^. 

In many-particle systems the Lyapunov exponent is 
also an indicator (order parameter) of phase transitions 
[^|-^ and may be related to transport coefficients . 

In practice, for a given system, the largest Lyapunov 
exponent must be obtained through numerical simula- 
tions, typically using the method developed by Benet- 
tin and collaborators @] (a proof that this method gives 



the Lyapunov exponent of Oseledec's theorem |Q| can be 
found in Ref. |jl^). For some special cases, e.g., hard- 
sphere systems, the theory of Lyapunov exponents is re- 
markably developed jl^]. This contrasts with the situa- 
tion in smooth Hamiltonian systems, where comprehen- 
sive analytical estimates are scarce. In 1984, when study- 
ing two-dimensional billiards, Benettin made the first 
step towards the construction of random matrix meth- 
ods for modeling universal features of Lyapunov spectra 
|l2| . These methods were further developed by other au- 
thors and applied to many-particle smooth Hamiltonians 
II- 

To our knowledge, the first theory for estimating the 
largest Lyapunov exponent of a specific Hamiltonian sys- 
tem, and its dependence on the system parameters, was 
formulated by Pettini and coworkers some years ago 
jl^ , ^^ . In this approach, the dynamics is geometrized 
by absorbing the force terms into a suitable metric, thus 
mapping the Hamiltonian problem onto a geodesic mo- 
tion on a curved manifold. After making the "quasi- 
isotropy" approximation, the Hamilton equations for the 
tangent vectors become decoupled. As a consequence, 
the initial system of 2N differential equations has been 
reduced to only two equations. While the original prob- 
lem was governed by the Hessian matrix of the potential, 
of size N X N, the new (reduced) one is controlled by 
the Laplacian of the potential, AV{t), a scalar function 
of time. Thereafter, AV{t) is treated as Gaussian white 
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noise and the 2x2 system of differential equations is solved 
using the methods developed by van Kampen and others 
p6t . See Ref. for a review. 

When applied to a Fermi-Pasta- Ulam chain the 
so-called "geometric method" was extremely successful in 
reproducing the largest Lyapunov exponent over the en- 
tire energy range [|l5|. However, in other cases the agree- 
ment is not so good. For instance, in a chain of rotators 
with first-neighbor (bounded) interactions, the method 
works well only in the low and high energy regimes, where 
the dynamics is weakly chaotic (integrable in the limits 
E ^ Q, oo). In the intermediate region of stronger chaos, 
the theory has to be amended to obtain a good agree- 
ment with simulations [p^ . This and other examples 
raise several questions concerning the domain of 
validity of the theory. What is the nature of the quasi- 
isotropic approximation? Or, what are the parameters 
that control the quality of the estimates of the theory? 
Is the geometric method perturbative? If so, what are 
the next leading corrections? 

In this paper we present an alternative theoretical ap- 
proach in which the validity domains of the successive 
approximations can be precisely delimited. The basic 
idea is to employ van Kampen's methods to solve 
the original system of 2N differential equations for the 
evolution of tangent vectors. By applying this scheme 
to a three dimensional dilute gas, Barnett et al. es- 
tablished a link between the Lyapunov exponent and the 
self-diffusion coefficient (see also jl^). We show that this 
approach can be extended to reach other systems, like 
Fermi-Pasta-Ulam chains and lattices of classical spins, 
either with short or long-range interactions. By doing 
so, we shall settle down a connection with the results of 
the geometric method and suggest some answers to the 
above-mentioned questions. 

The paper has been organized as follows. Section || 
presents the theory that leads to an estimate of the 
largest Lyapunov exponent. This is a perturbative the- 
ory, which rests on a cummulant expansion. We argue 
that the general (perturbative) solution can be obtained 



N 



by diagonalizing a small-dimension matrix. In Sect. Ill 
we analyze an approximation that re duce s the problem 
to diagonalizing a 3x3 matrix. Section IV discusses some 



examples which illustrate the working of the theory. The 
connection between our results and those obtained by 
the geometric method is discussed in Sect. 0. Finally, 
Sect. contains the concluding remarks. 



II. THEORY 



(1) 



where qi and pi, are conjugate position- momentum coor- 
dinates. Other Hamiltonians can be considered but they 
may require modifications of the theory. 

The Hamilton equations can be written in the compact 
form 



dx 



(2) 



where we have introduced the 2iV-dimensional column 
vector X, X = (gi, . . . ,qN,Pi, ■ ■ ■ ,Pn)'^, the superscript 
meaning "transposed", and the symplectic matrix J, 



1 

-1 



(3) 



with 1 the N X N identity matrix. Differentiating the 
Hamilton equations, one obtains the evolution equations 
for tangent vectors ^ = {Sqi, . . . , Sq^, Spi, 



,SpnV, 



(4) 



For a Hamiltonian of the special form (|l|), and setting 
m — 1, the operator A has the simple structure 



Mt) 
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-V{t) 



(5) 



Here V is the Hessian matrix of the potential V, namely 



dqidqj 



(6) 



Once initial conditions and have been specified, 
Eqs. (H) and (^) allow one to find the Lyapunov expo- 
nent A by calculating the limit jsj 



A - lim _ln|e(t;a;o,^o)|' 

t^oo Zt 



(7) 



We will assume that for any initial condition xq, the 
phase-space trajectory xit^xo) is ergodic on its energy 
shell. This implies that A is independent of initial condi- 
tions xo, which can then be chosen randomly according 
to the microcanonical distribution. There will also be 
no dependence on initial tangent vectors, because if is 
also chosen randomly, it will have a non-zero component 
along the most expanding direction. 

If the corrections to the exponential law in Eq. (^ go 
to zero fast enough as t — > oo, one can also write 



The theory we present in this section can, in princi- 
ple, be applied to any smooth Hamiltonian system. For 
simplicity, and for the sake of comparisons with the geo- 
metric method, we restrict ourselves to "natural" Hamil- 
tonians, 



(|C(t;xo,eo)|')(x 



„2\t 



(8) 



brackets meaning microcanonical averages over xq. We 
will prefer the estimate of Eq. (||) because the averaging 
procedure is crucial for finding an analytical expression 
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for the Lyapunov exponent. In case of doubt, the equahty 
of the exponents defined by Eqs. (0) and (||) can be tested 
numerically, e.g., using the data generated by Benettin's 
algorithm. 

By letting xq be a random variable, a stochastic pro- 
cess Y{t;xo) is defined, and Eq. (||) can be thought as 
a stochastic differential equation. However, the quantity 
we are interested in is the square of the norm of ^, which 
can be written as the trace of the "density matrix" 
Thus, we must focus on the equation for the evolution of 



A = i Re (L„ 



(15) 



^A^ EE 



(9) 



the rightmost identity defining the linear superoperator 
A. Except for the fact that we must deal now with a su- 
peroperator, Eq. (^) is not different from Eq. (||) and can 
be handled with the same techniques. For the purpose 
of the perturbative approximations that will follow, the 
operator A is split into two parts 



A = Ao + Ai(t) 



(10) 



where Aq corresponds to the evolution in the absence of 
interactions. In our case Aq and Ai are associated with 



1 




and Ai 











-V{t) 



(11) 



respectively. Whenever Ai{t) is small (in a sense that 
will be discussed below), it is possible to manipulate 
Eq. (^ to derive an explicit expression for the evolu- 
tion of the average of A clear exposition of this 
derivation, together with a very detailed discussion of its 
domain of validity has been given by van Kampen [[l6[ . 
We just outline the basic steps: (a) Rewrite Eq. in 
the interaction representation associated with Aq. (b) 
Write the propagator as a time ordered exponential, (c) 
Expand its average in a series of cummulants. (d) Go 
back to the original representation. The final result is: 



(12) 



where A is a time-independent superoperator given by 
the perturbative expansion: 

A = Ao + (Ai) 

dr (SAiit) e^^« SAi{t - r) e^^^M + • • • , (13) 

with 



In Eq. ([1^) we give explicitly only the first two cum- 
mulants, the dots stand for third cummulants and higher 
order ones. The perturbative parameter can be under- 
stood as the product of two quantities. The first one, 
let's call it a, characterizes the amplitude of the fluctua- 
tions of SAi{t). The second, Tc, is a typical (the largest 
relevant) correlation time of 6Ai{t). Thus, the second 
cummulant is of the order of ct^Tc, the third one is of 
the order of ct^t^, and so on. If all cummulants were 
summed up, Eq. (|l^) would be exact in the long-time 
regime Tc . 

From now on, we restrict our analysis to the propaga- 
tor A truncated at the second order, i.e., Eq. ( p^ ) without 
the dots. This approximation will be better the smaller 
CTTc. However, if Ai(i) is not far from a Gaussian pro- 
cess, the validity of the second order approximation may 
extend outside the perturbative region arc <C 1. In the 
exceptional case that Ai(t) is a Gaussian process, cum- 
mulants higher than the second one will be strictly zero 
and the truncation will introduce no error. 

To proceed further one needs the matrix of A in some 
basis. So, let us calculate AM, with M a symmetric ma- 
trix (it is easy to see that the truncation has not spoiled 
the symmetry of the density matrix). First notice that 
the exponentials of Aq represent no problem as they are 
finite polinomials, 

e^^oQ = [l + rAo]Q[l + TA;^], (16) 

for any matrix Q. Inserting this expression into (^|) we 
arrive at 

AM = (Ao + (Ai))M 

driSAiit) [6Ai{t-T)M + MSA^ {t-r)]] 

+ {---f, (17) 
where (• • ■)'^ means "the previous terms transposed", and 

6Ai (t-T) = [a + tAo] <5Ai {t - t) [a - tAo] . (18) 

Substituting (^Tj) into ( p7| ) we arrive at the final result of 
the second-order perturbative approach 

AM=( ;W 



, , ■ \ / {6YSV') , , . 



5Ai(t) = Ai(0- (Ai) 



(14) 



Let Linax be the eigenvalue of A which has the largest 
real part. We find that the largest Lyapunov exponent A 
is related to the real part of Lmax: 



dr 



0\ SY' 
SV \ <5V' 



— T — r 



(19) 
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To abbreviate the notation, we have written r" instead 
of r"l; 5Y and SV substitute SY{t) and SY{t - r), 
respectively. 



The largest Lyapunov exponent is buried into Eq. (19). 
To get an explicit expression one must diagonalize the 
matrix of A. The outcome will be A as a function of 
the first two cummulants of the stochastic process V(t), 
i.e., averages and (integrated) two-times correlation func- 
tions: 



dTT"{6V,j{Q)SVki{T)) , 71 = 0,1,2. (20) 



At first sight it may be thought that, as A is a superop- 
erator, the matrix one should diagonalize is of the order 
of N"^ X , then straightforward diagonalization would 
be out of the question for large N. Notwithstanding, A 
is an averaged object, and, as such, it possesses some 
symmetries which can be exploited to reduce the dimen- 
sionality of the problem to tractable levels, say N x N. 
This will be illustrated with an example in Sect. [V, So, 
if desired, the largest Lyapunov exponent could be found 
by numerical diagonalization, at least for systems with 
N « 1000 degrees of freedom (provided one can estimate 
the correlation functions). 

An alternative to exact diagonalization is approximate 
diagonalization, i.e., the diagonalization of the restriction 
of A to some small-dimension subspace. Notice that the 
problem we are dealing with is not very different from 
finding the ground state energy of a quantum system. 
However, in the quantum problem, the operator one must 
diagonalize is Hermitian, and it is well known that diago- 
nalization in a truncated basis produces an upper bound 
to the exact ground-state energy. It frequently happens 
that this bound is close to the exact result, even if the 
ground-state wavefunction is not. In spite of the operator 
A not being Hermitian (see Sect. [H) we still expect that 
diagonalization in a small basis will give a lower bound 
for the Lyapunov exponent. If the basis is suitably cho- 
sen, this estimate may be close to the result of the exact 
diagonalization. 

To proceed with the construction of a basis for A, we 
take advantage of the fact that A is independent of ^q, 
and simplify Eq. (|l^) further by averaging over an or- 
thonormal set of initial tangent vectors, obtaining 



(21) 



This second averaging allows us to consider, instead of A 
itself, the restriction of A to the subspace spanned by the 
matrices A'^1, A; = 0, 1, 2, . . .. A look at the first terms of 
these sequence gives a hint for constructing an appropri- 
ate basis. The first term is the identity, the second one 
is 



Al = 





1- (V) 



1-(V) 








t{SVSY') 



t{SVSV') (l-r2)((5V5V') 



(22) 



and so on. 



III. THE ISOTROPIC APPROXIMATION 

Typically, the diagonal elements of V(t) will be larger 
than the off-diagonal ones. This is evident in the case of 
translational invariancc, where one has the property 



(23) 



Introducing a matrix Y having all entries equal to one, 
i.e.. 



Yy=l, VZ,J. 

we can rewrite Eq. (^^ as 

YV = VY = . 



(24) 



(25) 



Then it is clear that Eq. (^3|) is also satisfied by (V), 
{SV5Y'), and by the higher moments of V that will ap- 
pear in the blocks of A*^!, for A: > 1. So, in a first 
(crude) approximation one may be tempted to discard 
the off-diagonal elements of the moments of V. If we also 
assume that all coordinates qi are statistically equivalent, 
and remind that the matrices A^^l are symmetric, we ar- 
rive at the simplest approximation for diagonalizing A. 
We call this approximation "isotropic" , and it consists in 
restricting A to the subspace spanned by the following 
three matrices: 



Ii 



a 





I, 




1 



1 

a 



(26) 



These matrices are mutually orthogonal with respect to 
the standard Euclidean scalar product, i.e.. 



Tr(I,lJ) cx 4 



(27) 



Then the matrix elements of A with respect to the basis 
{Ii, 12,13} are 



A 



_ Tr([AI,-]lf) 
Tr(Ijf) 



(28) 



Using Eq. ( [19D and skipping some simple algebra, we ar- 
rive at the 3x3 matrix 





A" = I 2a^TP 

-M + 2a2ri'^ 

with the definitions 



2 
-2a2ri'^ | , (29) 

1 -2aVi^) 
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M = ^Tr(V) , 
^2 = iTTiiSVf) , 

Jo 



(30) 
(31) 
(32) 



where we have introduced the normahzed correlation 
function /(r) 



fir) 



1 



Na- 



■Ti-{SV{0)SV{t)) 



Na 



1 ^ 

-2 E('^^'^(0)'^^»^-M) 



(33) 



It is evident from Eq. (g9|) that the operator A is not 
Hermitian. Normahzation of the basis Ij wih not make 
A^^ symmetric. 

In the isotropic approximation the Lyapunov exponent 
is expressed in terms of the set of four parameters fi and 
^2^ik+i) ^ fc = 0,1,2. The parameters fj, and a are, re- 
spectively, the mean and variance of the stochastic pro- 
cess V(i), and can in principle be obtained analytically 
by calculating the corresponding microcanonical aver- 
ages. (In practice the calculations can be done in the 
canonical ensemble, and then connected with the micro- 
canonical results by the formula of Lebowitz, Percus and 
Verlet |§.) 

The characteristic time rj^^ is naturally interpreted as 
the correlation time of the process V(t). Its calculation 
requires the knowledge of the autocorrelation functions 
of Vij{t), which are system dependent. Moreover, the 
correlation functions of a given system will in general de- 
pend on energy. However, if the functional form of the 
correlation function /(r) is known (or conjectured), its 
parameters can also be calculated as thermal averages. 
For instance, if /(r) is approximately Gaussian, 



fir) 



-7T 



(34) 



the expansion of (V(O)V(r)) around t — gives an ex- 
plicit formula for the correlation time, namely 



1 




1/2 



(35) 



In this case ri^-* and Tc^^ are trivially related to Tc^^: 



.(2) 



2 r 



-(3) ^ 1 



-(1) 



(36) 
(37) 



(k) 

A purely numerical calculation of Tc may be very dif- 
ficult, if not impossible, because correlation functions es- 
timated from finite-length time series usually fail to damp 



as expected plj. Perhaps, a more sensible approach to 

(k) 

the estimation of Tc should start with a numerical study 
of the correlation functions; then a functional form for 
/(r) could be proposed, based on the short-time behavior 
of the numerical correlation functions; finally, the param- 
eters defining /(t) would be calculated as suitable ther- 
mal averages. An alternative, more powerful approach, 
involves the use of "memory functions" ||2^: They are 
related in a one-to-one way to correlation functions and 
seem to be more amenable to simple approximations (see, 
e.g., Ref. H). 



IV. EXAMPLES 

In this section we analize the application of the per- 
turbative theory of Sect. || to some simple models. We 
remark that, in principle, the theory is expected to be 
successful only in regimes where the Lyapunov exponent 
is very small. All the systems considered below exhibit 
regimes with vanishingly small Lyapunov exponents. It 
is understood that our discussion will be restricted to 
such regimes. 



A. Mean field XY-Hamiltonian 

Let us begin by analizing one special case in which the 
isotropic approximation of Sect. [II is exact. Consider 
the one-dimensional Hamiltonian l24|-EQ| 



JV N 

2^ ' 2N 



cos 



-Oj)] 



(38) 



This is the so-called mean-field XF-Hamiltonian. It rep- 
resents a lattice of classical spins with infinite-range in- 
teractions. Each rotator is restricted to the unit circle 
and it is therefore described by an angle < 9i < 2tt and 
its conjugate angular momentum Li, with i = 1, . . . , N. 
At the critical energy Ec — 3N/ 4 there is a second-order 
phase transition separating a disordered regime {E > Ec) 
from an ordered one {E < Ec). 

In both limits E' 0, cxd the Lyapunov exponent goes 
to zero. For a fixed energy E > Ec, A also goes to zero 
when N —^ 00. This behaviour has also been observed in 
a metastable disordered phase with E < Ec- The pertur- 
bative approach should be a good approximation in these 
regimes. Moreover, we argue that the infinite-range in- 
teractions justify the isotropic approximation. 

All single-particle averages are equal, and, given that 
the forces are independent of the distances between spins, 
all two-particle averages must also be equal. So, one has 



Cl, 
C2, 



(39) 
(40) 
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Notice that translational symmetry, Eq. (]23|), implies 
that 



C2 



Cl 



N- 1 



(41) 



This is the reason why the isotropic approximation will 
work in this case, i.e., ofF-diagonal matrix elements are 
indeed smaller than diagonal ones But let us keep 
the discussion quantitave, and rewrite Eq. as 



(V) =cil + C2(Y-l) 



(42) 



with Y defined in Eq. (|2J) . Using the time-reversal sym- 
metry of the stochastic process V(t), one can also show 
that 



(5V(0)(5V(r)) ^c[l + 4(Y - a) 



(43) 



Then all blocks of Al [Eq. (p^j belong to the subspace 
spanned by 1 and Y. Taking into account Eq. ( ^ ) and 



Y^ iVY 



(44) 



we conclude that the blocks of all the sequence be- 
long to the subspace {1, Y}. Thus, the relevant subspace 
for diagonalizing A is six-dimcnsional. It is spanned by 
Ii,l2,l3 [Eq. M)] and Yi, Y2, Y3, with the definitions 



Yi 



Y 






Y 



Y3 = 



Y 
Y 



(45) 



However, it can be shown (see Appendix) that the largest 
eigenvalue of the corresponding 6x6 matrix coincides 
with that of the isotropic 3x3 matrix, up to corrections 
of order In this way we have proven the validity 

of the isotropic 3x3 approximation for one-dimensional 
systems with infinite-range forces. 



B. Dilute gases 



Consider now a one-dimensional gas with Hamiltonian 



1 " 



N 



J2 

2J = 1 



1j) : 



(46) 



where q and p are linear coordinates, and we assume that 
the potential 1/ is bounded. For large enough energies, 
this system is disordered and weakly chaotic. All parti- 
cles, and all pairs of particles, are statistically equivalent. 
Then this problem is formally equivalent to the infinite 
range Xy-Hamiltonian: The isotropic 3x3 approxima- 
tion becomes exact. 

Of course, the statistical equivalence also holds for a di- 
lute three-dimensional gas with short-range interactions. 
In this case, Barnett et al. have shown that the largest 
Lyapunov exponent is found by diagonalizing a 4x4 ma- 
trix M. 



C. aXy-Hamiltonian 

There are cases in which no strong reasons exist to be- 
lieve that the isotropic approximation will work satisfac- 
torily. Consider, for instance, the arbitrary-range analog 
of the XF-Hamiltonian [EtI,^: 



1 ^ 



i=l 



1 

2N 



N 

E 



1 — cos( 



(47) 



The parameter a sets the range of the interactions: a = 
recovers the mean-field case, a = 00 corresponds to 
first-neighbors couplings. The prefactor TV (a function 
of N and a) is included to make the system "pseudo- 
extensive" 29 1 . Periodic boundary conditions are as- 
sumed, and r.ij is the minimum between \i — j\ and 
N — \i — j\. For any value of a there exist (i) a low- 
energy regime of harmonic oscillators weakly coupled by 
non-linear forces and (ii) a high-energy disordered phase 
where the spins rotate almost freely. We expect our the- 
ory to produce good estimates for the largest Lyapunov 
exponent in both low- and high-energy regimes. 

If the forces are not infinite-range but just long-ranged, 
the isotropic approximation will still give good estimates 
in weakly chaotic regimes. Evidence supporting this 
statement can be found in Refs. [^ (geometric method) 
and [ |3l[ | (random matrix approach), where some kind of 
"isotropic" approximations were used to derive scaling 
laws for A in the high-energy regime, in good agreement 
with numerical simulations [ pTp^ ]. 

For a not too small, it may be necessary to improve the 
isotropic approximation by diagonalizing A in a larger 
basis. In this case, the statistical equivalence holds for 
all pairs of particles separated by the same distance. This 
means that the blocks of A''! are symmetric and cycli- 
cal, i.e., the matrix elements only depend on the distance 
Vij . A basis can be constructed starting from the NxN 
matrix S of a cyclical shift: 



(48) 



where it is understood that j + 1 must be taken modulo 
N. Then the set of symmetrical matrices 



< fc < N/2 , 



(49) 



is a basis for the blocks of A™!. A suitable basis for 
diagonalizing A is the set 











s, 











(50) 



with <i,j,k < N/2. The length of this basis is 3iV/2. 
Notwithstanding, we expect that a small subset of this 
basis will be enough to get a satisfactory convergence to 
the largest eigenvalue of A. Even in the worst case of no 
truncation at all, numerical diagonalization is possible 
for relatively large systems. 
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V. CONNECTION WITH THE GEOMETRIC 
METHOD 



In Sect. Ill we motivated the isotropic approximation 



by arguing that, in a first approach, one can neglect the 
ofF-diagonal matrix elements of the blocks of Then, 
in Sect. we proved that this approximation is indeed 
justified in various cases. Looking back to the results of 
Sects. H and [11, we realize that the isotropic approxima- 
tion is equivalent to postulating an "effective" system of 
equations. 





-K{t) 
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(51) 



where — (5qi,6pi) is the projection of the tangent 
vector ^ on the subspace of the i-th degree of freedom. 
The equations above represent the evolution of a typi- 
cal component of ^. In this sense they could alterna- 
tively be called "mean-field" or "single-particle" equa- 
tions. The scalar object K{t) is an effective random pro- 
cess which substitutes the Hessian V{t), and is in prin- 
ciple unknown. However, its first two cummulants can 
be identified in the following way. First solve Eq. (|5l| ) 
for the average of by using second order perturbation 
theory, as done in Sect. || [just change V hy K, and set 
A^=l in Eq. (|l9|)]. Notice that, as K{t) is a real num- 
ber, the blocks of the effective A are also real numbers, 
and the "isotropic approximation" is exact now. Then 
the matrix one must diagonalize to obtain the Lyapunov 
exponent is exactly that of Eq. (|2^) , provided one makes 
the identifications: 

{K) = iTr(V) , (52) 

{5K{Q)SK{t)) = lTr(<5V(0)5V(T)) , (53) 

with SK = K — (K). From this point of view, the 
isotropic basis {Ii, 12,13} is a single-particle basis. It 
is the most natural one in the sense that it treats all 
degrees of freedom on the same footing. 

The perturbative-isotropic approximation, as pre- 
sented above, is very similar to the geometric approach. 
In fact, in the geometric method an effective equa- 
tion like (^ij) is proposed, containing the unknown pro- 
cess K'it). Then, it is argued that the first two cummu- 
lants of K' are related to the Laplacian of the potential, 
/\V{t): 



(54) 



{5K'{Q)5K'{t)) = 1 {{5AVf) f5{r) , (55) 

where (5AV = AV — (AV), and f is the correlation 
time of the process K'{t), which is assumed to be delta- 
correlated. 



It is obvious that the averages of both processes K 
and K' coincide because TrV = AV. So, the differences 
between the geometric method and the perturbative- 
isotropic approach appear only in the fluctuations. The 
similarity between both theories could be enhanced by re- 
laxing the delta-correlation assumption of the geometric 
method, and substituting Eq. (ES) by 



1 



{SK'{Q)SK'{t)) = —{SAV{0)6AV{t)) 



(56) 



But even so, we have not been able to find any ana- 
lytical relationship between the correlation functions of 
K and K': In principle, the difference between both is 
non-negligible, and both effective theories will lead to dif- 
ferent estimates for the Lyapunov exponent. We expect 
that numerical simulations will decide which estimate is 
better. 

One comment about the correlation time f of Eq. ( [55| ) 
is in order: Geometric arguments lead to the estimate 



(57) 



with p. = {K') and 5-2 = {{dK')^). Some slightly different 
expressions have also been proposed ^,^,0. The crite- 
rion for testing the accuracy of theses estimates has been 
the agreement between the geometric estimate for A and 
numerical simulations, i.e., the goodness of fit (which is 
indeed excellent in some cases). To our knowledge, there 
is no independent test of the expression (|5^), or others, 
in the literature. Accordingly, a precise definition of f 
seems to be lacking. (Is f equal to the integral of the 
normalized autocorrelation function of AV{t)7) This is 
a point which affects the consistency of the geometric 
method. Unless a definition is given, to some extent, f 
will have the status of a fitting parameter. Comparisons 
of the geometric method with other theories will have to 
take this fact into account. 



VI. SUMMARY 

We showed that the evolution equation in tangent 
space can be thought of as a stochastic differential equa- 
tion with multiplicative noise. Then, an analytical es- 
timate for the largest Lyapunov exponent of a many- 
particle system in equilibrium was derived by using stan- 
dard perturbative techniques. Our analysis has been fo- 
cused on the second-order approximation. In this case 
the Lyapunov exponent can be obtained by diagonal- 
izing a matrix whose entries are calculated from the 
first two cummulants of the Hessian of the potential en- 
ergy, i.e., the averages (Vij) and the correlation func- 
tions {SVij {0)SVki{T)) . The dimension of this matrix is, 
in principle, of the order of NxN, but we have proposed 
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the conjecture, based on an analogy with the Hermitian 
problem, that diagonalization in a truncated basis may 
be enough to obtain satisfactory results. 

In the crudest approximation, which consists in choos- 



ing the three-dimensional isotropic basis of Eq. (26), the 
Lyapunov exponent is extracted from a 3x3 matrix. We 
argued that this "isotropic approximation" is equivalent 
to modeling the tangent dynamics of the many-particle 
system by an "effective" process K{t) for a single degree 
of freedom. In this way we established a connection with 
the so-called geometric method, the alternative effective 
theory for estimating the Lyapunov exponent. Both the- 
ories are very similar, but differ at the point of the defi- 
nition of the correlation function of K{t). The difference 
is non-trivial and is expected to lead to different predic- 
tions. 

In special cases, e.g., one-dimensional lattice-systems 
with infinite-range interactions, we have been able to 
prove that the isotropic approximation is exact. How- 
ever, in the general case, it may be necessary to consider 
larger bases. We have given examples where these bases 
are constructed by following the symmetries of the mo- 
ments of V(t). 

The theory we have presented is perturbative. Loosely 
speaking, we expect to obtain good estimates of Lya- 
punov exponents in weakly chaotic regimes. More quan- 
titatively, the domain of validity of the theory is con- 
trolled by the "Kubo number" arc, which quantifies the 
strength of the fluctuations <5V(i). For a given system, it 
is difficult to say a priori in which regimes the theory will 
work satisfactorily. This question and others, like the va- 
lidity of the isotropic approximation, and its comparison 
with the geometric method, will be decided with the aid 
of forthcomming numerical simulations. 
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The 6x6 matrix of A can be naturally split into four 
blocks of size 3x3. The block I-I has already been cal- 
culated [Eq. Let's now calculate the block I-Z, i.e.. 



^^/z^ Tr(AZ,)I, 



TrI? 



(A2) 



By setting V = in Eq. (|T^) we obtain the operator Aq 
(A in the absence of interactions). It has the following 
properties: 



AY, = AoYj , 
Tr ( AoY,) I, = Tr f AoL, ) I 



(A3) 
(A4) 



Using these two properties together with Eq. (Al) one 
arrives at 



KlZ _ J^II _ A // 



Analogously one obtains 

kZZ _ A// , }_KlI 

where the symbol ~ means that terms of relative size 
have been discarded. The 6x6 matrix reads 



(A5) 

(A6) 
(A7) 



(A8) 



with the definition Aq^ + A{^ = A-^^. Then it can be 
checked that the matrix above has three zero eigenval- 
ues while the remaining three are the eigenvalues of the 
matrix 



N+1 



II 



,11 



(A9) 



Thus the isotropic approximation is essentially exact for 
the infinite-range Xy-Hamiltonian. 



We are grateful to A. O. Caldeira, C. H. Lewenkopf, 
F. Nobre, A. M. Ozorio de Almeida, H. M. Pastawski, A. 
Robledo, A. M. C. de Souza, and C. Tsallis for fruitful 
discussions. We acknowledge Brazilian Agencies CNPq, 
FAFERJ and FRONEX for financial support. 



APPENDIX A: THE INFINITE-RANGE CASE 

We have seen in Sect. |^ that in the case of a 
one-dimensional system with infinite range interactions 
the subspace spanned by the matrices A*^! is six- 
dimensional. An ortoghonal basis for this subspace is 
the set {Ii, I2, 13, Zi, Z2, Z3}, where 
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